//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "ephi.hpp"

// see calculation noise on the axis of a ring of current/charge

#define RADIUS 0.1
#define WIDTH 800

int main (int argc, char *argv[])
{
  Statics statics;
  ElectroDynamics dyn(statics);

  ShapedCoilFactory(1).make (statics, vect3d(0, 0, 0), vect3d(1, 0, 0), RADIUS, 0.01, 100000, 1e-7);

  if (1)
    {
      Screen screen(WIDTH, WIDTH, 1e-2, 1e-2);
      Scene scene(dyn);
      screen.set_axis (Screen::PLANE_YZ, 100);
      scene.calc (screen);
      scene.set_coloring(32, true);
      scene.render_map (screen);
      screen.write ("test13_axis_bfield.bmp");
    }

  if (1)
    {
      Screen screen(WIDTH, WIDTH, 9e-13, 9e-13);
      Scene scene(dyn);
      screen.set_axis (Screen::PLANE_YZ, 0.000001);
      scene.calc (screen, Scene::CALC_EFIELD);
      scene.set_coloring(32, true);
      scene.render_map (screen);
      screen.write ("test13_axis_efield.bmp");
    }

  if (1)
    {
      prec_t w = 1;
      while ( w <= 100 )
        {
          vect3d tbf;
          dyn.getBField(vect3d(w, 0, 0), tbf);
          prec_t d = 1e-7;
          while ( d > 0 )
            {
              vect3d bf;
              dyn.getBField(vect3d(w, d, 0), bf);
              if ( fabs(1 - bf.length()/tbf.length()) > 0.01 )
                  break;
              d -= 1e-12;
            }
          if ( d < 0 )
              printf ("OK: %f\n", prec2double(w));
          else
              printf ("ERROR: %f, %e\n", prec2double(w), prec2double(d));
          w = w + 1;
        }
    }

  if (1)
    {
      Statics statics;
      ElectroDynamics dyn(statics);

      statics.addStaticElement (new StaticLine(100000, vect3d(-0.1, 0, 0), vect3d(0.1, 0, 0), 0.001));

      Screen screen(WIDTH, WIDTH, 1e-3, 1e-3);
      Scene scene(dyn);
      screen.set_axis (Screen::PLANE_YZ, 100);
      scene.calc (screen);
      scene.set_coloring(32, true);
      scene.render_map (screen);
      screen.write ("test13_line_bfield.bmp");
      scene.set_coloring(1);
      scene.render_lic (screen);
      screen.write ("test13_line_lic.bmp");
    }

  return 0;
}

